

def Ln(x):
  from scipy import special
  return(special.log(x))

def Wu(u):
  """
  Funcion que calcula Wu en base a desarrollo en series
  Parameters
  -----------------
  u : double
    Donde u es r^2*S/(4*T*t)
  """
  if u <= 0: 
    return ( [0, 0] )      
  eps = 0.0001
  su = - 0.577215665 - Ln(u)
  if u <= eps:
    return ( [su, u] )      
  tt = (-1)*u
  su = su - tt      
  if u < 1:
    n = 2
    while abs( tt ) >= eps:
      tt = (-1)*tt*u*(n-1)/(n*n)
      su = su - tt 
      n = n+1
      # print "n=%d,u=%f,tt=%f,su=%f" % (n, u, tt, su) 
    return ( [su, u] )  
  else:
    if u<=21:
      # Con N ~ 4*(u+1) la convergencia es bastante buena
      for cont in range(2, 4*(u+1)):
        tt = (-1)*tt*u*(cont-1)/( float(cont*cont) )
        su = su - tt      
        # print "cont=%d,u=%E,tt=%E,su=%E" % (cont, u, tt, su) 
      return ( [su, u] )
    else:
      # Arroja valores muy bajos y es demasiado costoso calcularlo
      # Ej. u = 18 Wu = 6.768184E-10
      return ( [eps*0.00001, u] )

def semillaT(tt, s, r,  Q,  Si = 0.1,  Ti = 100.0):
  print "tt=%s\ns=%s\nr=%f\nQ=%f" % (tt, s, r, Q)
  ui = []
  wui = []
  dd = []
  coef = Q / (4*3.1416*Ti)
  print coef
  cont = 0
  for t in tt:
    ui.append( r*r*Si/(4.0*t*Ti) )
    wui.append ( Wu( ui[cont] )[0] )
    print "ui=%f,wui=%f" % ( ui[cont], wui[cont] )
    ddd = s[cont] - coef * wui[cont]
    dd.append(  ddd )
    # print ui[cont], wui[cont], dd[cont]
    cont = cont + 1
  su = sum( dd )
  print "Suma incial = %f" % su
  if abs ( su ) <= 0.1:
    return ( Ti )
  if su >= 0:
    b = 5
  if su <= 0:
    b = 0.2
  b0 = b
  # su0 = su
  while True:    
    suant = su
    bant = b0
    su = 0
    if b>=1:
      for cont in range(0, len(tt)):
        su = su + ( s[cont] - Q * ( wui[cont]/b + 0.2 ) / ( 4 * 3.1416 * Ti ) )
    else:
      for cont in range(0, len(tt)): 
        su = su + ( s[cont] - Q * ( wui[cont] *b- 0.2 ) / ( 4 * 3.1416 * Ti ) )          
    print "\tb0 = %f, Suma = %f" % (b, su)    
    b0 = b
    if abs ( su ) <= 0.01:
      print "Transmisividad estimada=%f" % (b*Ti)
      return ( b*Ti )
    if su > 0:
      if b < 1:
        b =  b*1.1
      else:
        b = b*0.9
    if su < 0:
      if b < 1:
        b =  b*0.9
      else:
        b = b*1.1
      
      
if __name__ == '__main__':
  # Custodio Confinado
  r = 10.0
  ti= [1,3,6,10,20,40,70,100,200,400,700,1000,2000,3000]
  t = []
  s = [3.4,4.2,4.8,5.08,5.60,6.05,6.40,6.65,7.10,7.50,7.90,8.25,8.70,8.90]
  Q = 100.0/1000.0*60*60*24.0
  for tt in ti:
    t.append(tt/24.0/60.0)
  semillaT( t, s, r,  Q, Si = 0.3, Ti = 10.0 )
  # semillaT( t, s, r,  Q, Ti = 90.0 )
  # semillaT( t, s, r,  Q, Si = 0.05,  Ti = 700.0 )
